Defect self-propulsion in active nematic films with spatially varying activity

We study the dynamics of topological defects in active nematic films with spatially varying activity and consider two set-ups: (i) a constant activity gradient and (ii) a sharp jump in activity. A constant gradient of extensile (contractile) activity endows the comet-like +1/2 defect with a finite vorticity that drives the defect to align its nose in the direction of decreasing (increasing) gradient. A constant gradient does not, however, affect the known self-propulsion of the +1/2 defect and has no effect on the −1/2 that remains a non-motile particle. A sharp jump in activity acts like a wall that traps the defects, affecting the translational and rotational motion of both charges. The +1/2 defect slows down as it approaches the interface and the net vorticity tends to reorient the defect polarization so that it becomes perpendicular to the interface. The −1/2 defect acquires a self-propulsion towards the activity interface, while the vorticity-induced active torque tends to align the defect to a preferred orientation. This effective attraction of the negative defects to the wall is consistent with the observation of an accumulation of negative topological charge at both active/passive interfaces and physical boundaries.


Introduction
Active nematics are collections of elongated apolar particles that consume energy from their surroundings to generate dipolar forces that drive self-sustained flows [1]. Much progress in understanding the rich dynamics of these active liquid crystals has been achieved through a minimal hydrodynamic theory that couples orientational order and flow and captures the behaviour of biological systems from subcellular to multicellular scales [2].
What distinguishes the hydrodynamics of active nematics from that of their passive counterparts is the presence of an active stress generated by active processes, which sets up spontaneous spatiotemporally chaotic flows [7]. The active stress is given by s a ij ¼ aQ ij with Q the nematic order parameter and α a scalar activity parameter that encapsulates the biochemical processes that generate active forces [10][11][12][13]. It can have either sign: α > 0 corresponds to a system of 'pullers' generating contractile stresses on their surroundings, whereas α < 0 reflects a system of 'pushers' and their induced extensile active stresses. With increasing activity, active flows are induced spontaneously and create large distortions of the nematic order, including the formation of pairs of topological defects that sustain active turbulence [3,14].
The lowest energy topological defects in active nematic films have half-integer charge, corresponding to the comet-shaped +1/2 with polarity e + defined by a head-tail arrow, and the −1/2 which has threefold symmetry (see figure 1). These defects disrupt the nematic order locally and induce longrange distortions in the orientation field, generating active stresses, which in turn lead to spontaneous active flows surrounding the defects [12,15,16]. There is a net active flow through the core of the polar +1/2 defect which makes it intrinsically motile and is referred to as the defect self-propulsion. For an isolated +1/2 defect, the self-propulsion velocity aligns with the polarity vector and, depending on the contractile/extensile properties of the active nematic, the defect moves in/opposite to the direction of its polarization. The −1/2 defect does not create any net flow at the defect position and, thus, is not self-propelled in systems with uniform activity. The motion of defects in the presence of spatially inhomogeneous activity is far less understood and explored [17].
There are several approaches to realize experimentally systems with spatially dependent activity. In [18], a varying substrate topography is used to control the frictional damping in a film of a microtubulekinesin suspension. This results in spatial variations of the concentration of active agents, thus indirectly in the local activity. In [19], a similar effect is achieved by manipulating light-sensitive myosin motors that activate the microtubules. Both studies find that the −1/2 defects localize near the interface separating the region of higher activity from that of lower activity. In [19], it was reported that the +1/2 defects are deflected by the active/passive interface. Analytical work based on a hydrodynamic theory of the defect gas has predicted that a passive/active interface can be used to separate positive and negative topological charge [17]. Numerical studies of how the defect dynamics is affected by the spatially dependent activity show that the polarity of the +1/2 defect tends to align parallel to the activity gradient [19][20][21][22][23], and that the confinement and motion of defects can be manipulated by varying the steepness of the activity gradients [24,25]. A recent numerical study also shows that the formation of defect dipoles can be controlled by imprinting special geometries into the activity profile [22].
In this paper, we provide a theoretical study of how spatially varying activity affects the selfpropulsion and reorientation of isolated topological defects. We consider the representative basic set-ups where the spatial activity profile is given either by a constant activity gradient (linear profile) or a sharp interface separating two regions of constant bulk activity (Heaviside function profile). For constant activity gradients, the +1/2 defect rotates due to a vorticity-induced active torque acting on the defect polarization until the defect aligns parallel to the activity gradient and moves in the direction of lower magnitude of activity. Thus the defect slows down. We show analytically that the vorticity at the +1/2 defect core is proportional to the hydrodynamic dissipation length ' d ¼ ffiffiffiffiffiffiffiffi ffi h=G p , which measures the strength of viscous dissipation η relative to friction G. Numerical simulations of royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 the active flow generated in a disc of radius R show that the vorticity depends on the system size for small R, and crosses over to the analytically predicted value for large systems. The vorticity field induced by an activity gradient parallel to the +1/2 defect's polarization has a quadruple structure with regions of alternating vorticity. This is confirmed by numerical simulations for a disc geometry where four vortices are formed around the +1/2 defect. By contrast, the vorticity induced by constant activity gradients at a −1/2 defect has an eightfold symmetry that leads to eight vortices with alternating circulation in a finite domain. We also calculate the net translational self-propulsion and reorientation that both ±1/2 defects acquire near a sharp active/passive interface. The +1/2 defect slows down as it moves towards the interface, and the vorticity-induced torque tends to reorient it such that its polarization becomes normal to the interface regardless of the sign of activity. The −1/2 defect also acquires a preferred orientation at the interface, and those that approach the interface with this stable orientation are then attracted by it. The structure of the paper is as follows. We start in §2 by introducing the minimal hydrodynamic model active nematic films on a substrate. In §3, we derive and discuss the self-propulsion and spontaneous rotation of +1/2 defects in the presence of a constant activity gradient. Section 4 focuses on the analytical derivation of the self-propulsion and vorticity of ±1/2 defects close to a sharp active/passive interface. Summary and concluding remarks are presented in §5.

Hydrodynamics of active nematics with spatially varying activity
We consider the familiar hydrodynamic model of a 2D active nematics that couples the flow velocity uðrÞ to the nematic order parameter Q ij ¼ Sðn inj À 1 2 d ij Þ, where S quantifies the degree of order and nðrÞ ¼ ðcos uðrÞ, sin uðrÞÞ is the orientational director field with head-tail symmetry. In the simplest formulation, the Q-tensor is a minimizer of the de Gennes-Landau free energy [7] with isotropic elastic constant K > 0 and g the strength of the local ordering potential. The uniform nematic ordered state corresponds to S 0 = 2. The flow field satisfies a Stokes equation that balances forces on a fluid element, given by [7] (G À hr 2 )u ¼ r Á ½aðrÞQðrÞ À rpðrÞ, r Á u ¼ 0, ð2:2Þ where G is a friction coefficient per unit area, η is the shear viscosity and α is the activity coefficient, with dimensions of stress. For simplicity, we neglect the contributions from the passive stresses and flow alignment to focus, instead, on the active flows generated by an isolated ±1/2 in the presence of nonhomogeneous activity aðrÞ.
We rescale the Stokes equation in units of the nematic relaxation time τ = γ/g (where γ is the nematic rotational friction) [16,26] and the coherence length j ¼ ffiffiffiffiffiffiffiffi ffi K=g p . Different dynamical regimes are then controlled by one dimensionless , and the rescaled activity aðrÞ ! aðrÞg=ðGKÞ. The dimensionless form of the Stokes equation reads as where the active force field induced by an isolated ±1/2 defect is given by The first contribution is an interfacial force F I + originating from activity gradients. The second term is a bulk force F B + due to nematic distortions. The defect self-propulsion velocity v + is defined as the net active flow through the defect core, and thus can be computed from the active flow velocity u obtained from the solution of equation (2.3) evaluated at the origin [26]. Due to viscosity, it depends nonlinearly and non-locally on the force field through the integral solution of equation (2.3) given by where K 0 (r) is the zeroth order Bessel function which is the Green's function of equation (2.3) without the incompressibility constraint. We distinguish the interfacial contributions v I + from the bulk contributions v B + . The incompressibility constraint gives rise to pressure gradients which may affect the defect royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 3 kinematics. The pressure field is the solution of the corresponding Poisson's equation The net vorticity at the defect core is also obtained from measuring the vorticity of the flow field induced by the defect distortion, given by v ¼ @ x u y À @ y u x ¼ Àr ? Á u. Using equation (2.5) and evaluating it at the defect position r 0 ¼ 0, we obtain an expression for the defect vorticity Vorticity is also written as sums of interfacial v I + and bulk v B + contributions which depend on the defect polarization e + and are computed analytically in the next sections.
For isolated ±1/2 point-like defects, we can parametrize the Q-tensor order parameter in the quasistatic phase approximation as [26] Q + xx ðrÞ ¼ cosð+fðrÞ þ 2u 0 Þ and Q + xy ðrÞ ¼ sinð+fðrÞ þ 2u 0 Þ, ð2:8Þ where fðrÞ ¼ arctanðy=xÞ is the singular part of the nematic orientation due to a ±1/2 defect located at the origin, and θ 0 is the slowly varying part of the background orientation of the nematic director. The +1/2 defect has a well-defined polarization which is determined by the background nematic orientation θ 0 as For the −1/2 defect, we can also introduce a polarization vector determined by θ 0 and aligning with one of the principal axes of the threefold symmetry [26] e À ¼ cos 2u 0 3 , sin 2u 0 3 : ð2:10Þ Both nematic defects and their respective polarizations are illustrated in figure 1. It can be shown that a net vorticity at the defect core induces an active torque that tends to rotate the defect polarization. This follows straightforwardly from taking the time derivative of the polarization in equations (2.9) and (2.10), and using the evolution of the Q-tensor [17,26] to account for the change in the background nematic field θ 0 due to vorticity as ∂ t θ 0 ≈ ω/2. Thus, the evolution of the defect polarization controlled by vorticity is where the defect charge is q = ±1/2 and e ? ¼ ½e y , À e x represents the 90°clockwise rotation of the polarization vector. For motile defects, there are additional torques due to defect interactions, the elastic stiffness K or the coupling to the flow alignment [17,26]. Here, we focus on the active torque induced by a non-zero vorticity which emerges from spatially varying activity alone. In the subsequent sections, we investigate how this active torque reorients the defect polarization relative to activity gradients for two set-ups: (i) a constant activity gradient and (ii) an interface with a sharp jump in activity.

Constant activity gradient
We first study the kinematics of an isolated defect in a region where the activity gradient is locally constant. Without loss of generality, we consider an activity gradient in the x-direction such that the activity has the linear profile α(r) = α 0 + α g x. The defect orientation is arbitrary and controlled by the background nematic orientation θ 0 . We demonstrate that a constant gradient α g does not modify the defect self-propulsion velocity as compared with what was obtained for uniform bulk activity α 0 . An activity gradient across the texture of a +1/2 defect generates, however, a flow that may yield a finite vorticity at the defect core, which tends to align the defect polarization according to equation (2.11) in the direction of the gradient. The −1/2 defect remains stationary both in its motion and orientation.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 The interfacial active force given in equation (2.4) arising from a constant activity gradient α g is wherer is the radial unit vector andr ? ¼ ðŷ, ÀxÞ. This expression corresponds to the Helmholtz decomposition of F I þ into a curl-free part ( r) and a divergence-free part ( r ? ). These two contributions are plotted in figure 2. The divergence-free part gives a net vorticity at the defect core which tends to rotate its polarization until it aligns with the activity gradient. This is most easily demonstrated in the friction-dominated limit where the active flow velocity is Gu ¼ F I þ À rp. The incompressibility constraint thereby removes the curl-free contribution through the contribution of the interfacial pressure which is radially symmetric and given by p I þ ðrÞ ¼ a g cosð2u 0 Þðr À LÞ, ð3:2Þ making the interface flow purely rotational. Here the constant L is a length comparable with the system size which controls the divergent terms. More generally, to incorporate viscous dissipation we need to evaluate the integral expression for the defect velocity given in equation (2.5). In an infinite system, the symmetry of the integrand leads to no contribution to the defect speed from the interfacial active force, thus v I þ ¼ 0. This contribution may become finite in non-radially symmetric bounded domains. The contribution from the bulk active force in equation (2.4) reduces to where F 0 þ ðrÞ is the known active force corresponding to a constant activity α 0 which leads to a constant self-propulsion velocity [16,26]. Since the contribution due to activity gradient α g is antisymmetric royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 around the defect position, its integral according to equation (2.5) vanishes. Therefore, there is no contribution from activity gradients to the defect self-propulsion. This is not changed when we add the gradient of the bulk pressure which is given as Here p 0 þ is the pressure for the constant activity term α 0 [16]. Activity gradients induce, however, a vortical flow that is finite at the defect core, resulting in an angular velocity of the +1/2 defect, given by where the first term in the bracket originates from the interfacial active force and the second is due to the bulk force. The integral can be carried out in polar coordinates, with the result v þ ¼ 3pa g 4z sinð2u 0 Þ: ð3:6Þ We can rewrite this equivalently in physical units as to highlight that the defect angular velocity scales linearly with the hydrodynamic dissipation length l d , similar to the self-propulsion speed of a defect in a constant activity [16]. The effect of this vorticity is to align the polarization so that it is pointing opposite to the activity gradient. This is consistent with recent numerical results, where defects align normal on soft interfaces separating extensile and contractile regions [23].
To test the validity of these analytical predictions for a bounded system, we have solved numerically the Stokes flow from equation (2.3) in a disc of radius R. Equation (2.2) is solved with non-slip boundary conditions using the finite-element package FEniCS [27,28]. The active stress is computed from the analytical form of the Q tensor corresponding to a single point defect in a uniform background nematics.
In figure 3, we show that the defect angular velocity is proportional to R for radii smaller than l d , and crosses over to the asymptotic value for an infinite system given by equation (3.7) at large R. We have also computed the vorticity field for α 0 = 0 and different defect orientations relative to the activity gradient, as shown in figure 4. When the defect polarization is parallel to the activity gradient (θ 0 = 0), we observe a quadruple structure of the vortical flow. This is consistent with the analytical prediction in the frictiondominated limit, where the vorticity field away from the defect is determined by the activity gradient α g as (for α 0 = 0) v þ ðr, fÞ ¼ a g sinð2fÞ 2Gr cosð2u 0 Þ þ a g Gr sinð2u 0 Þ(1 þ cos 2 ðfÞ), ð3:8Þ where r and ϕ are the polar coordinates centred at the defect position. By contrast, when the defect polarization is normal to the activity gradient (θ 0 = π/2), we obtain a single vortex centred at the core of the defect.   royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229

−1/2 defect
In the friction-dominated limit and for α 0 = 0, we can evaluate the vorticity field, and show that a constant activity gradient α g induces eight counter-rotating vortices, with a vortical flow given by v À ðr, fÞ ¼ À a g cosð2u 0 Þ 2Gr ð3 sinð4fÞ À sinð2fÞÞ þ a g sinð2u 0 Þ 2Gr (3 cosð4fÞ À cosð2fÞ): ð3:11Þ The same structure is observed in bounded domains where the vorticity forms vortices of alternating circulation, as shown in figure 5 for a disc geometry.

Activity jump at an interface
We now consider an activity profile corresponding to a sharp interface separating a region of high activity α 0 from a region of low activity α 1 . Isolated ±1/2 defects are situated at a distance x v from the interface in the region of high activity, α 0 , as illustrated in figure 6. The activity profile across this interface is given by the Heaviside step function corresponding to a singular activity gradient ∂ x α = −Δαδ(x − x v ) with Δα = α 0 − α 1 the interfacial jump in activity. An active/passive interface corresponds to α 1 = 0 and Δα = α 0 . In this case, we find that the selfpropulsion of the +1/2 defect is reduced as the defect approaches the interface. The vorticity-induced active torque tends to reorient the ±1/2 defects moving toward the interface to preferred orientations that depend on extensile/contractile activity. The −1/2 defect that already has the selected orientation is attracted to the wall, while that with different polarizations might be repelled.  x v x v e + a 0 a 1 a 0 a 1 Figure 6. Set-up of (a) +1/2 defect and (b) −1/2 defect at a sharp interface separating a region with higher activity α 0 from that with lower activity α 1 .

+1/2 defect
The active force field induced by a +1/2 defect located at a distance x v from a sharp interface is given by Inserting these expressions in equation (2.5), we obtain the contributions to the self-propulsion velocity from interfacial and bulk active forces as The first term in the bulk contribution is the well-known constant self-propulsion velocity from a constant activity α 0 [16]. The second term is the additional drift due to the activity jump Δα and depends on the distance x v from the interface. As we will see below, this contribution suppresses the defect self-propulsion near the interface.
If we now specialize to the case of an active/passive interface, i.e. Δα = α 0 . In dimensional units, the self-propulsion velocity of the +1/2 defect is then given by The function f þ v ðx v Þ is plotted in figure 7a. We note that the self-propulsion speed vanishes as the defect hits the interface x v = 0. In other words, the defect slows down as it approaches the interface, and eventually remains at rest at the interface. We note that equation (4.5) is obtained by incorporating the incompressibility constraint only in the v 0 þ term. Additional pressure gradients may arise due to activity jump. These are, however, difficult to obtain analytically and are not included in this study.
We now compute the vorticity at the defect position to investigate how its contribution to the active torque tends to reorient the defect as it approaches the interface to a stable orientation. From equation  Figure 7. Plot of (a) f þ v and (b) f þ v as functions of the distance x v of the +1/2 defect from the interface. Note that f þ v diverges at x v = 0 due to the bulk terms.
royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 (2.7), we obtain the following expressions for interfacial and bulk contributions where the bulk vorticity diverges at x v = 0. The total defect angular velocity is given by the sum of these two contributions evaluated at the defect core. In dimensional units, it is given by where the wall-dependence function f þ v ðx v Þ is plotted in figure 7b. Hence, near an active/passive interface, the vorticity-induced rotation is at a rate _ It is clear from figure 7b that reorientation only occurs within a distance of order ℓ d from the wall. As the +1/2 defect approaches the wall, f þ v increases and eventually diverges at x v → 0. This means that the defect tends to reorient its polarization until sin(2θ 0 ) = 0. From the stability criterion that (dω + /dθ 0 ) < 0, this corresponds to the stable orientation 2θ 0 = π for α 0 < 0 (extensile) and 2θ 0 = 0 for α 0 > 0 (contractile). In both cases, the defect polarization is normal to the interface e þ ¼ ½+1, 0 and points away from the interface for extensile systems and into the interface for contractile systems, respectively. Numerical simulations [19] report that +1/2 defects tend to reorient and drift parallel to the boundary when the angle between the interface and the incoming velocity is below a critical value that depends on activity. Above this critical angle, i.e more head-on collisions, the defect hits the wall and tunnels through it. This effect is probably coming from the additional contributions to the active torque that are not considered here, namely the interactions between defects, deformations in the nematic order parameter due to the wall and the coupling to flow alignment. It is likely that these terms are important close to the interface, both for determining the defect orientation and the tunnelling effect observed both experimentally and numerically [19].

−1/2 defect
The components of the interfacial active force due to a −1/2 defect at a distance x v from the activity jump are given by The corresponding bulk active force is 1 r 3 ½ðy 2 À x 2 Þ cos 2u 0 À 2xy sin 2u 0 ð 4:12Þ and F B yÀ ¼ ða 0 À DaHðx À x v ÞÞ 1 r 3 ½ðy 2 À x 2 Þ sin 2u 0 þ 2xy cos 2u 0 : ð4:13Þ Using these expressions, and neglecting the contribution from the pressure gradient, the net drift velocity of the defect can be written as royalsocietypublishing.org/journal/rsos R. Soc. Open Sci. 10: 221229 wheren À ¼ cosð2u 0 Þx þ sinð2u 0 Þŷ. The function f À v ðx v Þ describes the dependence on the distance x v to the interface and is given by dx dyK 0 ðrÞ x 2 À y 2 r 3 : It has been evaluated numerically and is plotted in figure 8a. The −1/2 defect acquires a finite selfpropulsion close to the wall in a region of thickness of order ℓ d near the activity jump. Its motion is either towards or away from the boundary, depending on the defect's orientation and the sign of the activity. To see how the −1/2 reorients as it approaches the interface, we evaluate the flow vorticity at the defect core as a function of the wall distance. Again, there are contributions to the vorticity from both flows driven by interfacial and bulk forces, given by ð4:16Þ Note that the bulk term diverges when x v → 0. The total angular velocity of the −1/2 defect can then be written as The function f v ðx v Þ has been calculated numerically and is shown in figure 8b. The dependence on the wall distance x v changes sign near the wall, indicating that the vorticity tends to rotate the defect to a preferred orientation at the wall. The preferred orientation is determined by the stationary condition sin2θ 0 = 0, and the stability criterion (dω − /dθ 0 ) < 0, which implies that θ 0 = 0 for α 0 < 0 and 2θ 0 = π for α 0 > 0. In other words for extensile activity the stable orientation of a −1/2 defect at a sharp active/passive interface corresponds to a polarization e − = [1,0]. Therefore, as a result of both their self-induced translational and rotational motion, in an extensile system −1/2 defects are attracted to a sharp active/passive interface and orient themselves with one of the three axis normal to the interface. This is consistent with the accumulation of negative topological charge observed in experiments at active/passive interfaces [18,19] and near physical walls [29], as well as in simulations [24,25].

Conclusion
Activity gradients or sharp jumps can guide the motion and orientation of nematic defects. In a constant activity gradient, +1/2 defects acquire an angular velocity that may rotate their orientation such that the defect polarization aligns parallel to the activity gradient. The defects then self-propel in the direction of the gradient, always moving towards regions of lower magnitude of activity, where it is less motile. Thus, we expect that activity gradients will introduce more circular motion in the trajectories of the +1/2 defects. By contrast, a constant activity gradient yields no net vorticity or active force at the core of the −1/2 defect, which remains stationary. We find that the self-propulsion velocity of +1/2 defects moving towards a sharp active/passive interface is also reduced, and that the defect will eventually stagnate at the wall. By contrast, −1/2 defects acquire a finite propulsion speed in the interfacial region and can overcome the positive defects, explaining the observation of negative charge accumulation in experiments and simulations [18,19,24,25,29]. We also predict that the active torque acting on a +1/2 defect that reaches the interface tends to reorient it toward a preferred polarization that is perpendicular to the interface and points away/toward it depending on extensile/contractile activity. The vorticity-induced active torque also acts on the orientation of a −1/2 defect migrating toward interface, by rotating the defect until it reaches the stable orientation which minimizes the net vorticity at the defect position. We show that a −1/2 defect with a stable orientation gets attracted to a sharp interface. This stable orientation is selected by the sign of activity, i.e whether the system is contractile or extensile. Tunnelling across the interface observed numerically may be due to soft interfaces where the activity gradients are not sufficiently steep, as well as due to defect interactions and other hydrodynamic effects. Here, we have neglected additional contributions of pressure gradients induced by activity gradients, as well as elastic stresses, flow alignment, nematic distortions due to the active/passive interface and defect interactions, which may change qualitatively the defect dynamics.
Our results offer a simple understanding of the dynamics of nematic defects in the presence of spatially varying activity. They can provide the starting point for designing structures capable of controlling defect dynamics and associated active flows.